function res = residual_for_labor(Par,phi,vz,agg,vloglabor,zguess)
% Description:
%   Returns the residuals of the first order conditions of firm's profit
%   maximzation problems at given level of labor allocation, taking
%   productivity and aggregate variables as given
% Output:
%   res - 4 by 1 vector of residuals
% Input:
%   Par - parameter structure
%   phi - scaler, parameter phi for the Home firm
%   vz - 2 by 1 vector of productivity ([zHome;zForeign])
%   agg - 1 by 3 vector of equilibrium values of aggregate variables
%   vloglabor - 4 by 1 vector of log of labor input
%   zguess - scaler, initial guess of cutoff productivity in log
% Note: Labor allocation in log form (vloglabor) to accommodate searches of
%   negative numbers by fsolve
% Note: Productivity and labor input CANNOT be matrices; Only vectors (for
%   individual sectors) are allowed for this function
    
    % Load and rename paramters
    kappa = Par.kappa;
    sigmazH = Par.sigmazH;
    vlabor = exp(vloglabor);
    zbarH = log(vz(1,1))-sigmazH^2/2;
    zF = vz(2,1);
    wF = agg(3);
    % Turn off display of message by fsolve
    options = optimoptions(@fsolve,'Display','off');
    % Solve for cutoff productivity level ztildeH
    objfun = @(logztildeH) residual_for_ztilde(Par,phi,vz,agg,vlabor,logztildeH);
    z_sol = fsolve(objfun,zguess,options);
    ztildeH = exp(z_sol);
    % Define a function of marginal revenue product of labor to be used in
    % the integrals of the FOCs
    fmrpl = @(zH) fun_mrpl(Par,agg,[zH;zF],vlabor);
    % Integrals for Home firm's first order conditions
    intH1 = @(zH) (fmrpl(zH) - fmrpl(ztildeH) + (1-phi).*(fmrpl((1-kappa)*ztildeH)-(1-kappa))).*lognpdf(zH,zbarH,sigmazH);
    intH2 = @(zH) (fmrpl((1-kappa)*zH)-(1-kappa)).*lognpdf(zH,zbarH,sigmazH);
    % Assign residuals from the Home first order conditions
    % Only 1st and 2nd row are used; 3rd and 4th row to be removed later
    res = integral(intH1,ztildeH,Inf,'ArrayValued',1) + (1-phi).*integral(intH2,0,ztildeH,'ArrayValued',1);
    % Integrals for Foreign firm's first order conditions
    intF1 = @(zH) (fmrpl(zH)-wF).*lognpdf(zH,zbarH,sigmazH);
    intF2 = @(zH) (fmrpl((1-kappa)*zH)-wF).*lognpdf(zH,zbarH,sigmazH);
    resF = integral(intF1,ztildeH,Inf,'ArrayValued',1) + integral(intF2,0,ztildeH,'ArrayValued',1);
    % Assign residuals from the Foreign first order conditions
    res(3:4) = resF(3:4);
end